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Abstract 

A wide range of organisms use circadian clocks to keep internal sense 
of daily time and regulate their behavior accordingly. Most of these clocks 
use intracellular genetic networks based on positive and negative regulatory 
elements. The integration of these "circuits" at the cellular level imposes 
strong constraints on their functioning and design. Here we study a recently 
proposed model [N. Barkai and S. Leibler, Nature, 403:267-268, 2000] that 
incorporates just the essential elements found experimentally. We show that 
this type of oscillator is driven mainly by two elements: the concentration 
of a repressor protein and the dynamics of an activator protein forming an 
inactive complex with the repressor. Thus the clock does not need to rely on 
mRNA dynamics to oscillate, which makes it especially resistant to fluctua- 
tions. Oscillations can be present even when the time average of the number 
of mRNA molecules goes below one. Under some conditions, this oscilla- 
tor is not only resistant to but paradoxically also enhanced by the intrinsic 
biochemical noise. 
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The environment changes in a highly periodic manner. There are, among 
others, daily cycles of light and dark as well as annual cycles of changing 
climates and physical conditions. Such environmental periodicity may create 
the necessity for organisms to develop internal time-keeping mechanisms to 
accurately anticipate these external changes and modify their state accord- 
ingly In particular, a wide range of organisms, as diverse as cyanobacteria 
and mammals, have evolved circadian rhythms — biological clocks with a pe- 
riod of about twenty four hours that evoke and regulate physiological and 
biochemical changes to best suit different times of the day. 

Recent findings show that the molecular mechanisms upon which these 
clocks rely share many common features among species p|. The main charac- 
teristic is the presence of intracellular transcription regulation networks with 
a set of clock elements that give rise to stable oscillations in gene expres- 
sion. A positive element activates genes coupled to the circadian clock. It 
simultaneously promotes the expression of a negative element, which in turn 
represses the positive element. The cycle completes itself upon degradation 
of the negative element and re-expression of the positive element. 

A crucial feature of circadian clocks is the ability to maintain a constant 
period over a wide range of internal and external fluctuations Such ro- 
bustness ensures that the clock runs accurately and triggers the expression 
of clock-dependent genes at the appropriate time of the day. For instance, 
fluctuations in temperature affect chemical reaction rates and may perturb 
oscillatory behavior. Another source of fluctuations may be the presence of 
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internal noise due to the stochastic nature of chemical reactions 0. Low 
numbers of molecules may be responsible for random fluctuations that can 
destabilize the oscillatory behavior of the biochemical network ||^. Yet, cir- 
cadian clocks maintain a fairly constant period amidst such fluctuations. 

Description of the model 

To study possible strategies, or principles, that biological systems may use 
to minimize the effect of stochastic noise on circadian clocks, we examine a 
minimal model based on the common positive and negative control elements 
found experimentally [|]. This model is described in Figure |l|. It involves two 
genes, an activator A and a repressor R, which are transcribed into mRNA 
and subsequently translated into protein. The activator A binds to the A 
and R promoters, which increases their transcription rate. Thus, A acts as 
the positive element in transcription, whereas R acts as the negative element 
by sequestering the activator. 

The deterministic dynamics of the model is given by the set of reaction 
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rate equations 



dDA/dt = OaD'^^-^aDaA 
dDn/dt = eRD'j,-^nDRA 
dD'Jdt = ^aDaA-OaD'j, 

dnydt = jrDrA-OrD'r 

dMA/dt = a'AD'j^ + aADA-huMA 
dA/dt = (3AMA + eAD'^ + eRD'ji 

-A{^aDa + IrDr + 7c^ + 5 a) 
dMR/dt = a'RD'ji + aRDR-6MnMR 
dR/dt = I3rMr-^cAR + 5aC -5rR 
dC/dt = ^cAR-5aC, 

where the variables and constants are as described in the caption for Figure |l|. 
This simple model is not intended to reproduce the particular details of each 
organism but to grasp the properties that the core design confers. As in any 
general model, the parameters of the values we use are typical ones. For 
instance, the rates for bimolecular reactions are all in the range of diffusion 
limited reactions. 

The preceding equations would be strictly valid in a well-stirred macro- 
scopic reactor. At the cellular level, a more realistic approach has to con- 
sider the intrinsic stochasticity of chemical reactions ||5|. This can be done 
by transforming the reaction rates into probability transition rates and con- 
centrations into number of molecules. One then obtains the so-called master 
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equation which gives the time evolution of the probability of having a given 
number of molecules. There is no general procedure to solve this type of 
equation analytically, but it is the starting point to simulate the stochastic 
behavior of the system. The basic idea behind such simulations is to perform 
a random walk through the possible states of the system, which are defined 
by the numbers of molecules of the different reacting species. Starting from 
a state with given numbers of molecules, the probability of jumping to other 
state with different numbers of molecules (i.e. the probability for an elemen- 
tary reaction to happen) can be computed from the master equation. One 
can pick up a state and the jumping time according to that probability dis- 
tribution, then consider the resulting state as a new initial state, and repeat 
this procedure until some final state or time is reached. In this way, the num- 
bers of molecules change in time with the statistical properties given by the 
master equation. There are several algorithms to implement this. The main 
difference among them is the specific way in which they compute the prob- 
abilities and select the states. For chemical reactions with few components, 
it is customary to use the so-called Gillespie algorithm P|. 

This intrinsic probabilistic behavior in the evolution of the number of 
molecules gives rise to fluctuations that are usually referred to as noise. In 
general, the term noise is used for any disturbance interfering with a sig- 
nal or with the operation of system. In the case of chemical reactions, the 
signal would be the average production of the reacting species whereas the 
disturbance would arise as a consequence of the fluctuations around that av- 
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erage value. We use term noise rather than fluctuations to emphasize the 
disturbing efl'ect that these fluctuations can have. Thus, although related, 
both terms do not mean the same. For instance, there can be large fluctua- 
tions in some molecular species but, if their characteristic time is very short 
compared to those of other processes that take place, they would introduce 
little noise. 

In Figure H we compare the results of the stochastic and deterministic 
approaches. We show the levels of A protein and R protein over time for 
the set of parameter values and initial conditions given in the caption of Fig- 
ure |l[ The deterministic results were obtained from numerical integration of 
Eqs. (|l|), whereas the stochastic results were obtained by computer simulation 
using the Gillespie algorithm. The main difference between the determinis- 
tic and stochastic time courses is the presence of random fluctuations in the 
latter. In the deterministic model every circadian cycle is identical to the 
previous one. The stochastic model shows some variability in the numbers 
of molecules and the period length, corresponding to the intrinsic fluctua- 
tions of the biochemical network. For these values of the parameters, both 
stochastic and deterministic approaches give similar results. We have also 
used different initial conditions and in all the cases we have observed that 
the behavior of the long term solution is the same. 

Model simplification 

To gain further insight into the essential elements that are responsible for the 
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oscillations, we will simplify as much as possible the deterministic rate equa- 
tions. By making various quasi-steady state assumptions ||^, it is possible to 
simplify the set of Eqs. (|T]) into a two variable system with the repressor R 
and the complex C as the two slow variables: 



dt e^ + ^^A{R) ^ 

— = ^cAim-^AC (2) 



where 

A{R) = \{a'AP{R) - Kd) + yia'MR) - KdY + ^aAp{R)K, 

with p(-R) = Pa/ ^Ma^IcR + ^a) and = Oa/ia- Notice that the nonlin- 
earity in the equations enters trough the quasi-equilibrium value of A, A{R), 
which is a function of R. The main idea behind these approximations is that 
there are fast and slow variables. Fast variables are assumed to be at an 
effective equilibrium whereas slow variables are responsible for the dynam- 
ics of the system. Thus, given the set of Eqs. (|^), we assume that all the 
derivatives except dR/dt and dC/dt are zero. 

In Figure || we show that for the values of the parameters we use the 
numerical solutions for the trajectories of the two-variable system [Eqs. @] 
agree closely with the solutions of the full system [Eqs. ([![)], except for quan- 
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titative differences in the peak levels and times at the beginning of each 
cycle. These differences arise because the time scale separation between fast 
and slow variables is not sufficiently large for quasi-steady state assumptions 
to be exact. These results indicate, nevertheless, that the dynamics of the 
system is mainly determined by two component concentrations: those of the 
complex and the repressor. The other components are driven mainly by these 
two elements and their effects enter the system through effective parameters. 
It is worth emphasizing that the reduced two-variable model is aimed just 
to offer insights into the qualitative behavior of the system and to show how 
the properties that one observes in the full system are already present in a 
simple two variable model. Thus, whenever we present simulation results for 
the deterministic system, except if otherwise stated, we are referring to the 
full system. Regarding the validity of the two-variable model, it is a good 
approximation when the dynamics of mRNA and the activator is faster than 
that of the complex and repressor. For instance, it will remain valid if 6^ is 
decreased or some of 6a, 5a/a? ^Mr, Oa, and 9r are increased with respect to 
the parameters of the caption for Figure |l|. 

Limit cycle oscillations and stability analysis 

The existence of oscillations in the two- variable system can be inferred from 
application of the Poincare-Bendixson theorem f^. This theorem states that 
a two-dimensional system of the type we are considering exhibits limit cycles 
if it is confined in a closed bounded region that does not contain any stable 
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fixed points. The trajectory of the system is confined since the number of 
molecules cannot reach infinite values. The fixed points and their stability 
can be determined by following a standard linear stability analysis. There is 
a single fixed point for positive concentrations. In our case, the signs of the 
real parts of the eigenvalues of the matrix describing the linearized dynamics 
around this point are given by 



dA{R) 




OR 





bcA{R) + 5r + 5a] , (3) 



where r is the trace of the matrix. All the quantities are evaluated at the 
fixed point {Ro,Co). When r is positive, the real part of the eigenvalues is 
positive, the fixed point is unstable, and there is a limit cycle in the system, 
which gives sufficient conditions for the existence of oscillations. Evaluation 
of Eq. (H) shows that r is indeed positive for the set of parameters we are 
using. The domain in which r is positive is rather broad. For instance, the 
function r remains positive when 9a and 9r are multiplied by a factor K with 
0.024 < K < 10.7; when all transcription rates (a and a') are multiplied by K 
with K > 0.08; or when protein and mRNA degradation rates are multiplied 
by K with 0.0009 < < 3.5. When r is negative, the fixed point is stable 
and, in principle, no conclusion about the existence of limit cycles can be 
drawn. In this case, the presence of oscillations could also depend on the 
initial conditions. For the full model, the ranges of parameters that give rise 
to oscillations are not exactly the same but remain very close to the previous 
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ones. 

The mechanism responsible for oscillations is illustrated in Figure]^ through 
the phase portrait of the two-variable model. Starting with low numbers of 
initial molecules near the origin of the phase plane, the trajectory rapidly 
shoots upwards along the R = nullcline (the dot over a variable means its 
time derivative). Here, the high levels of A, present due to auto-activated 
transcription, rapidly induce the formation of the complex C. Reaching the 
maximum of the nullcline, the trajectory 'falls off' the edge and moves rapidly 
diagonally right and downwards, corresponding to a drop in C and a rise in 
R. The trajectory curves around the R = nullcline and hits the (7 = 
nullcline, where it slowly returns to the left and approaches the fixed point 
{Rq, Co). When approaching the fixed point, R decreases sharply, taking the 
trajectory past the fixed point and sending it back upwards to initiate a new 
cycle. 

The trajectory in Figure | comprises a fast phase corresponding to the 
rapid production of C and R, and a slow phase corresponding to the slow 
degradation of R. These two distinct phases are characteristic of excitable 
systems, the classic example of which is the Fitz Hugh-Nagumo model for 
action potential transmission in neurons [l^. The fast and slow leg cor- 
respond to the excitable and refractory phase of the system, respectively. 
Thus, the system oscillates as it avoids the fixed point {Ro,Co) and hits the 
R = nullcline on the left to begin the excitable phase of a new cycle. 

As we have already pointed out, the deterministic analysis can be useful 
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to grasp the main properties of the system under certain conditions. Unfortu- 
nately such conditions are not known a priori without a stochastic analysis. 
Surprisingly enough, we have found that parameter values that give rise to 
a stable steady state in the deterministic limit continue to produce reliable 
oscillations in the stochastic shown in Figure |. Therefore, the pres- 

ence of noise not only changes the behavior of the system by adding more 
disorder, but can also lead to marked qualitative differences. 

How can the system continue to produce oscillations even when deter- 
ministic rate equations predict a stable steady state? The system always 
evolves towards a stable fixed point, as sketched in Figure ||. However, a 
perturbation of sufficient magnitude near the fixed point, e.g. as illustrated 
by the dotted arrow in Figure || (notice that the figure has not been drawn to 
scale and that the size of the arrow is not representative of the actual size of 
the perturbation), may initiate a new cycle. For low numbers of molecules, 
the intrinsic fiuctuations of chemical reactions can be large enough to contin- 
ually send the system into the fast phase after each cycle and thus produce 
sustained oscillations. In the deterministic limit (or close to it), there are 
no perturbations (or the perturbations are too small) to initiate a new cycle 
and the trajectory stays close to the fixed point. In this case, the system 
performs better if enough noise is present in the system. This situation is 
analogous to that observed in the Fitz Hugh-Nagumo model, where an op- 



timal amount of noise maximizes the reliability of the oscillations ||1 1|| . It 
is important to realize that the effects that noise may have on non-linear 
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systems can be difficult to predict and rather paradoxical ||T^]. Therefore, 
the smaller number of molecules does not necessarily imply more irregular 
behavior of the system, as one might intuitively assume |^, |T^. 



Significant parameters and noise resistance 

The mechanism responsible for oscillations involves only two variables. This 
has some important consequences for the functioning of the clock. If we 
consider the deterministic limit, a two-dimensional dynamical system of this 
type either oscillates regularly or does not oscillate at all. In two dimensions, 
since trajectories cannot cross, fixed points and periodic orbits are the only 
possible attractors. Other more complicated behaviors such as chaos or quasi 
periodicity are not allowed ||^. On the other hand, the intrinsic stochastic 
fluctuations of the remaining variables are effectively averaged and do not 
significantly affect the performance of the system. 

For instance, one variable that usually plays a prominent role in many 



circadian rhythm models is the number of mRNA molecules In our 

case, however, mRNA does not enter directly into the dynamics. It is just an 
intermediate step in the production of the proteins. Thus, if protein produc- 
tion remains unaltered the system will oscillate regardless of the number of 
mRNA molecules involved. This can be achieved, for instance, by increasing 
simultaneously the degradation rates of mRNA and the translation rate of 
the proteins. In the deterministic limit of the two-variable model this has 
no effect on protein dynamics at all. In the stochastic case the effects are 
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negligible. Figure |^ shows the time evolution of repressor mRNA and protein 
levels in the system for (3 a-, Pr, ^Ma multiplied by a factor of hun- 

dred. The system essentially alternates between having zero and one mRNA 
molecule in the cycle, and the proteins continue to exhibit remarkably good 
oscillations. 

There are also parameters that do affect the properties of the oscillations. 
In the deterministic limit, oscillations are always regular (provided that the 
two- variable model is a good approximation). When fluctuations are taken 
into account, the reliability of the oscillations may depend on those param- 
eters. One such parameter, as we have seen in the previous section, is the 
repressor degradation rate 5r. This parameter affects the period of the os- 
cillations (compare e.g. Figures ^ and |5|) and also the stability of the fixed 
point. One can infer from Eq. (|^) that for high or low values of the fixed 
point becomes stable. In such cases, the deterministic system stops oscillat- 
ing but this does not need to be so for the stochastic one, which may continue 
producing reliable oscillations. 

Notice that the positive feedback is a key element in the clock dynamics. 
Its most obvious use is the generation of the instability that gives rise to 
oscillations. But it has other not so obvious role that is closely related to the 
resistance to noise. In general, gene regulation is a slow process (with typical 
characteristic times of about one hour) and, as a such, is prone to be affected 
by fluctuations. This problem gets even worse if the dynamics relies in several 
coupled transcriptional feedbacks since the effects of the fluctuations are then 
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amplified. The positive feedback gives a fast transcriptional switch, allowing 
to move fast from low to high transcription rates. In this way, the time in 
which the system is prone to fluctuations is greatly reduced. 

Conclusions 

The presence of noise in transcriptional and enzymatic networks is a funda- 
mental consequence of the stochastic nature of biochemical reactions. The 
ability to function effectively and consistently amidst such random fluctu- 
ations is a major issue in gene expression and network behavior. In this 
paper we have studied how different factors affect a simple model for circa- 
dian rhythms that exhibits noise resistance. We found that the oscillations 
in this model are driven mainly by two components: a repressor protein and 
an activator-repressor complex. This fact is responsible for the reliability of 
the oscillations. First, a two-dimensional dynamical system of this kind has 
a very simple behavior: in the deterministic limit, it either oscillates or goes 
to a steady state. Second, noise and perturbations in the other variables 
affect the system only slightly since they do not enter directly the dynam- 
ics. Finally, resistance to noise is achieved as the number of molecules of 
any of the two key components reaches small values only for short periods of 
time or when they are not driving the dynamics of the system. In this way, 
even though some molecular species may be present in very low numbers, the 
intrinsic stochasticity of biochemical reactions can be bypassed. 

It is important to emphasize that organisms have evolved networks to 
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function in the extremely noisy cellular environment. Suitable network de- 
signs, as those that are now emerging from the experimental data 15 



can 



confer resistance against such noise. In addition, some of these networks may 
not only be resistant to but could also be taking advantage of the cellular 
noise to perform their functions under conditions in which it would not be 
possible by deterministic means. 
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Figure Captions 

FIGURE 1: Biochemical network of the circadian oscillator model. 
and Da denote the number of activator genes with and without A bound 
to its promoter respectively; similarly, and Dr refer to the repressor 
promoter; Ma and Mr denote mRNA of A and R; A and R correspond to 
the activator and repressor proteins; and C to the the inactivated complex 
formed by A and R. The constants a and a' denote the basal and activated 
rates of transcription, /3 the rates of translation, S the rates of spontaneous 
degradation, 7 the rates of binding of A to other components, and 9 the 
rates of unbinding of A from those components. Except if otherwise stated, 
along this paper we have assumed the following values for the reaction rates: 
aA = 50 hr-\ a'^ = 500 hr-\ = 0.01 hr"\ = 50 hr"\ Pa = 50 hr"\ 
Pr = 5 hr-\ Sma = 10 hr-\ Smr = 0.5 hr-\ Sa - 1 hr~\ Sr = 0.2 hr-\ 
7A = 1 mor^ hr~\ 7ii = 1 moP^ hr'^, -fc ^ 2 moP^ hr"\ 9a = 50 hr~\ 
9r = 100 hr~^, where mol means number of molecules. The initial conditions 
are Da = Dr = 1 mol, D'^ = D'j^ ^ Ma = Mr = A ^ R ^ C ^ 0, which 
require that the cell has a single copy of the activator and repressor genes: 
Da + D'j^ = 1 mol and Dr + D'^ = 1 mol. The cellular volume is assumed to 
be the unity so that concentrations and number of molecules are equivalent. 
Notice that we assume that the complex breaks into R due to the degradation 
of A and therefore the parameter Sa appears twice in the model. 

FIGURE 2: Oscillations in repressor and activator protein numbers ob- 
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tained from numerical simulations of the deterministic (a, 6) and stochastic 
(c,(i) descriptions of the model. 

FIGURE 3: Time evolution of the quantities R (continuous line line) and 
C (dashed line) for the system reduced to two variables (a) by various quasi- 
steady state assumptions and for the complete system (6). 

FIGURE 4: Phase portrait of the two variable oscillator [Eqs. (||)] for the 
parameter values given in the caption for Figure |l] (the drawing is not to 
scale). The thick line illustrates the trajectory of system. {Rq.Cq) is the 
fixed point of the system, and R = dR/dt = and C = dC/dt = are the 
R and C nullclines respectively. The solid arrows give the orientation of the 
direction field on the nullclines. 

FIGURE 5: Time evolution of R for the deterministic [Eqs. (P) (a) and 
stochastic (b) versions of the model. The values of the parameters are as in 
the caption of Figure |l| except that now we set 5ji = 0.05 hr~^. For these 
parameter values r < so that the fixed point is stable. 

FIGURE 6: Phase portrait as in Figure | but for a situation in which the 
system falls into the stable fixed point {Ro,Co). The dotted arrow to the 
left of the fixed point illustrates a perturbation that would initiate a single 
sweep of the (former) oscillatory trajectory. 

FIGURE 7: Stochastic time evolution of the number of activator (a) and 
repressor (c) molecules; and activator (6) and repressor (d) mRNA molecules. 

19 



The values of the parameters are as in the caption for Figure [l| but now with 
Pa = 5000 hr-\ Pr = 500 hr"\ Sma = 1000 hr-\ and Smb. = 50 hr'K 
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